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ZPEMP:  A  ZERO-D  Z-PINCH  IMPLOSION  CODE 


I.  INTRODUCTION 

"The  discharge  of  electrical  energy  through  a  cylindrical  gas  puff  is  a 
well  tested  means  o£  creating  an  intense  burst  of  high  energy  radiation. 
Typical  output  from  an  imploding  neon  Z-pinch  is  a  few  kJ  of  x-rays  when  the 
peak  current  is  about  1  MA  on  the  GAMBLE  II  pulsed-power  device  at  the  Naval 
Research  Laboratory.  The  experimental  program  to  produce  a  bright  and 
consistent  Plasma  Radiation  Source  (PRS)  on  GAMBLE  II  and  on  other  facilities 
has  evolved  over  the  past  several  years..  This  program  rapidly  changes  as  new 

cdhcepts  in  plasma  load  design  are  continually  proposed.  To  predict  the 

utility  of,  analyze,  and  revise  such  new  ideas  inexpensively,  it  is  useful  to 
model  the  proposed  implosion  in  a  simulation.  A  zero-dimensional  radiation- 
hydrodynamics  code  provides  both  a  succinct  and  rapid  resource  to  address  the 
above  issues. 

In  a  zero-D  model  of  an  axisymmetric  implosion,  the  entire  plasma  is 

pictured  to  implode  as  a  unit,  either  in  an  annular  shell  or  a  filled 

cylinder.  The  velocity  of  the  imploding  shell  or  cylinder  is  determined  by 

the  volume  integral  of  the  momentum  or  kinetic  energy  equation;  the  internal 

energy  by  the  volume  integral  of  the  plasma  energy  equation.  Several  zero-D 

codes  have  been  used  in  the  past  to  predict  and  study  Z-pinch  dynamics.  The 

2  3  4 

code  SIMPLODE  was  used  to  study  neon  ’  and  sodium  implosions,  while  a 
similar  code  was  employed  by  McDonald  and  Ottinger^  to  also  study  a  neon  PRS. 

3 

In  the  SIMPLODE  code  the  mass  density,  velocity,  and  current  density  are 
taken  as  uniform  across  the  annular  shell.  The  waveform  of  the  driving 

current  waveform  is  a  prescribed  function  of  time,  not  determined  from  a 
circuit  equation.  The  radiative  cooling  is  based  on  a  detailed  non-LTE 
collisional-radiative  equilibrium  (CRE)  model  for  the  ionization  dynamics. 
The  escape  probability  approximation  is  used  for  the  radiative  transfer  of  the 
lines.  The  momentum  equation  determines  the  velocity  of  the  imploding  shell. 
Recent  investigations  with  the  SIMPLODE  code  indicate  a  serious 

Manuscript  approved  March  6.  1989. 
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drawback, however:  total  energy  is  not  always  conserved.  In  the  McDonald  and 
Ottinger  code^  the  mass  density  in  the  plasma  shell  is  uniform,  while  the 
velocity  profile  is  determined  from  the  continuity  equation.  The  velocities 
of  the  inner  and  outer  shell  radii  are  equal.  To  maintain  total  energy 
conservation,  the  outer  shell  velocity  is  determined  by  the  kinetic  energy 
equation  rather  than  the  momentum  equation.  The  current  density  is  uniform 
but  limited  to  a  small  skin  region  in  the  outer  part  cf  the  plasma  shell. 
Although  the  extent  of  this  region  of  resistive  heating  evolves  in  time,  is 
does  so  according  to  the  local  plasma  velocity.  Hence  the  same  mass  element 
is  always  carrying  the  current  and  no  diffusion  of  the  magnetic  field  is 
possible.  The  temperature  of  the  skin  is  allowed  to  be  different  from  that 
of  the  core  region.  The  ionization  state  and  radiative  emissivity  are 
determined  from  fits  to  a  sample  CRE  calculation.  Opacity  effects  on  the 
radiative  output  were  scaled  with  the  load  mass  to  match  experimental  data. 

The  present  ZPIMP  code  is  a  new  zero-D  model  which  adopts  the  best 
aspects  of  the  previous  codes,  namely  the  circuit  equation  and  CRE  ionization, 
while  improving  on  the  dynamics.  For  the  hydrodynamic  equations,  this  code  is 
structured  much  like  a  typical  multi-zone  1-D  simulation  code.  In  §11. A  to 
§11. C  the  spatial  and  temporal  finite  difference  forms  of  the  hydrodynamic 
equations  are  derived  in  detail.  These  results  are  reduced  to  three  zones  in 
§11. D,  with  most  of  the  mass  contained  in  the  central  zone.  This  zone  is 
characterized  by  an  averaged  density,  ion  temperature,  and  electron 
temperature.  The  finite  difference  approach  to  the  hydrodynamics  allows  one 
to  calculate  separate  velocities  for  the  inner  and  outer  radii  of  the  shell. 
The  steps  in  §11. D  also  point  out  the  limitations  of  a  zero-D  code, 
particularly  at  the  bounce  time.  For  the  evolution  of  the  magnetic  field,  the 
induction  equation  is  analytically  solved  (§11. E)  for  the  magnetic  flux  in  a 
moving  medium.  The  solution  assumes  a  small  gradient  length  scale  for  the 
field  and  uses  a  planar  geometry  for  the  region  of  magnetic  diffusion.  As  the 
plasma  heats  up  and  the  resistivity  decreases,  the  form  of  the  solution 
displays  the  onset  of  flux  freezing.  A  simple  circuit  model  is  employed  in 
§11. F  to  self  consistently  determine  the  driving  current  from  a  supplied  open 
circuit  voltage.  Section  II. G  addresses  the  problem  of  total  plasma  plus 
field  energy  conservation.  The  ionization  dynamics  is  handled  by  a  non-LTE 
CRE  model  (§III),  and  the  probabili ty-of-escape  method  for  radiative  transfer 
is  used  to  handle  opacity  and  trapping  effects  in  the  lines.  At  the  start-up 
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of  the  simulation,  where  the  temperatures  are  below  20  eV,  the  CRE  model  is 
replaced  by  an  average-atom  model  as  described  in  §IV.  These  last  two 
sections  are  brief  since  the  essential  features  have  been  described  in 
previous  papers. 

As  stated  in  the  opening  paragraph,  the  primary  motivation  for  a  zero-D 
simulation  code  is  the  rapid  investigation  of  new  concepts.  For  instance, 
recent  interest  has  turned  to  imploding  aluminum  vapor  puffs  on  the  GAMBLE  II 
device.  In  §V  we  apply  the  ZPIMP  code  to  model  such  implosions.  The 
implosion  radius,  maximum  electron  temperature,  maximum  charge  state,  and 
maximum  continuum  and  line  radiation  is  investigated.  The  initial  geometry  of 
the  vapor  is  varied  from  annular  puffs  to  filled  cylinders.  Detailed  studies 
of  various  mass  loadings  and  geometric  configurations  to  determine  the  optimum 
radiative  yields,  emitted  spectra  from  the  L-  and  K-shells,  and  pulse 
durations  will  be  presented  in  a  following  report. 

A  summary  of  the  essential  points  of  the  report  is  contained  in  §VI. 


II.  DYNAMICS 

A.  Magnetohydrodynamic  Equations  in  Cylindrical  Symmetry. 

The  structure  of  the  simulation  code  ZPIMP  is  modelled  after  the 
standard  form  of  large  scale  radiation-hydrodynamic  codes.  The  ionization 
dynamics  and  radiation  transport  is  time  split  from  the  magneto-hydrodynamics, 
and  communication  between  these  two  physical  aspects  is  handled  through  a 
local  radiative  heating  or  cooling  term.  The  ionization  and  radiation  aspect 
will  be  discussed  in  §§  III  and  IV.  In  the  present  section  we  consider  the 
magnetohydrodynamics  of  the  Z-pinch  and  derive  the  appropriate  equations  for 
numerical  simulations  by  starting  from  the  general  plasma-fluid  equations. 
Let  p.=m.n.  ( p  =m  n  )  be  the  ion  (free  electron)  mass  density:  eZ.  the  mean 
charge  state  of  the  ions;  p^^  (electron)  pressure;  and  (e^) 

the  ion  (electron)  internal  energy.  We  will  assume  cylindrical  symmetry  (no 
•♦•-dependence)  and  axial  invariance  (no  Z-dependence)  so  that  all  plasma  and 
field  variables  depend  only  on  the  time  t  and  the  radial  coordinate  r. 
Furthermore,  we  assume  for  the  mass  weighted  velocity  v  =  and  for  the 

magnetic  field  B  =  last  assumption  together  with  Ampere'  ■  law 
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implies  J  =  current  density  in  the  absence  of  displacement 
currents.  Thus  the  electrons  have  the  same  velocity  as  the  ions  in  the  radial 
direction,  vis.  v^^,  but  unlike  the  ions,  there  must  be  a  non-zero  electron 
velocity  in  the  axial  direction. 

We  now  write  down  the  plasma  and  field  equations  subject  to  the  above 
geometry.  The  mass  continuity  equation  is 


at 


(1) 


The  electron  continuity  equation  is  replaced  by  the  requirement  of  charge 
neutrality,  n  =  n.Z. . 

The  total  momentum  equation  is  found  by  adding  the  electron  and  ion 
momentum  equations  and  neglecting  the  electron  inertial  terms: 


*  r  lF<'‘i''R>  »  lF<Pi*V“vis>  ■ 


- 

C  Z  9 


(2) 


In  this  equation  is  an  artificial  viscosity  term  used  to  numerically 

smooth  the  conversion  of  kinetic  energy  to  thermal  energy  at  the  bounce. 
Section  II. C  on  finite  differencing  discusses  the  appropriate  form  for  Qyjg* 
The  equation  for  the  ion  internal  energy  is 


*  r  lF‘”l=i’'R>  *  <Pi*0vis>r  If<"r>  *  F  ”) 


In  the  present  application  we  will  be  considering  gaseous  vapors  so  that  the 

strong  coupling  effects  of  high  density  plasmas  can  be  ignored.  For  an  ideal 

perfect  ion  gas  e.  =  (3/2)k_T./m. .  The  thermal  heat  flux  is 
1  011 


(A) 


where  is  the  thermal  conductivity  given  by  Braginskii^  or  Epperlein  and 
Haines  .  Since  we  are  carrying  separate  ion  and  electron  temperatures  the 
energy  transfer  from  the  electrons  to  the  ions  is 


m  n 

Q.  =  3  —  —  k-(T.-T  ) 

le  m.  T.  B  1  e 

1  le 


(5) 
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where  t.  is  the  ion-electron  collision  time, 
le 

The  equation  for  the  electron  internal  energy  is 
*  r  r  +  ?  !?<■'%, R> 

-  -  »ie*  "l-’z  ♦  k  ■’z  lF<Ve>  •  <^> 

e 

The  specific  electron  internal  energy  is  composed  of  two  parts,  the  first 
component  is  due  to  the  free  electrons  and  the  second  is  the  electronic  energy 
of  excitation  and  ionization:  e  =  (3/2)k„T  /m  +  (p./p  )£  .  The  radiative 

heating  rate,  or  cooling  rate  if  negative,  is  contained  in  A.  The  quantities 
E  and  A  are  calculated  prior  to  each  hydrodynamic  step  and  kept  constant 

€X 

during  that  step.  The  electron  heat  flux  is 


"..R  ■  -  ^  .  (7) 

e 

The  thermal  transport  coeffcients  K, ^  p.,  and  a,  are  again  given  by 
Braginskii”,  and  we  have  expressed  the  resistivity  as  =  Oj^/e‘'ng  .  Also, 
sign(B^)  =  B^/|B^1. 

In  place  of  the  electron  momentum  equation  a  generalized  Ohm's  law  for 
the  plasma  is  derived  by  dropping  the  electron  inertial  terms  and  solving  for 
the  electric  field.  The  Z-component  is  given  by 


e 

To  complete  the  set  of  equations  we  have  Faraday's  law  for  the  evolution 
of  the  magnetic  field. 


3B.  3E„ 

— i  =  r 


ar 


and  Ampere's  law  without  the  displacement  current  for  the  current  density. 
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For  later  reference  ve  present  the  equations  for  the  plasma  kinetic 
energy,  plasma  total  energy,  and  field  energy.  Multiplying  the  total  momentum 
equation  by  and  using  the  mass  continuity  equation  leads  to  the  plasma 
kinetic  energy  equation; 


3,1 

3F<2'>i''R  >  • 


1  3  r  fi  2 

r  ^  r''Rl2Pi''R 


Pi" 


Pe" 


Q-l  -  (p-+P  +  0-)~  f— (rv_) 
visjj  'i  ^e  vis^  r  3r'  R^ 


=  -  rh%  ■  <") 

Adding  to  this  result  both  internal  energy  equations  gives  the  total  plasma 
energy  equation; 


3t 


p.  e.  + 
k  1  1 


p  e  + 
e  e 


1 

2"i 


) 


1  3  f 

F  ^hi- 


1  9  f  f  1  2 

*■  ~  p.e.+  p  e  +  ;rp.v„  + 

r  3r  L  Rill  ^e  e  2^i  R 

f  1 

A  +  J.,  -  -VoB. 

Z I  c  R  ik  1  Z  en 


Pi"  Pe"  ^vis]] 


(12) 


Poynting's  theorem  for  the  total  magnetic  field  energy  is  readily  derived  from 
Faraday's  law  by  multiplying  by  B^/An: 


1  <^±] 

3t  I8nj 


7^  -  |-(rE„B.) 
4/1  r  3r  Z  ♦ 


■^z^z  ■ 


(13) 


The  above  equations  contain  all  the  transport  terms  as  given  by 
Braginskii  subject  to  the  conditions  ol  cylindrical  symmetry  and  translational 
invariance  along  the  axial  direction.  Following  previous  studies  we  will 
lienceforth  drop  the  thermal  frictional  component  of  the  electric  field  (the 
term  in  P^)  in  comparison  with  the  resistive  component  (the  term  in  hj^)  in 
eqns  (6),  (8),  and  (12).  For  ~  I/2iiR6,  where  I  is  the  current,  R  is  the 
plasma  radius  and  S  is  an  effective  skin  depth,  the  neglect  of  the  frictional 
component  is  valid  as  long  as 


I  » 


620  3/2 

Z.ln(A)  lp°j‘e 


amps . 
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In  this  condition  temperatures  are  in  eV,  distances  in  cm,  and  is  a 

factor  of  order  unity  depending  on  (eB^/m^c)Tjg  and  Z^. 


B.  Integral  Form  of  the  Equations. 

In  a  zero-D,  cylindrically  symmetric,  radiation-hydrodynamic  simulation 

code  the  basic  assumption  is  that  most  or  all  of  the  plasma  is  confined  to  an 

annulus  or  filled  cylinder.  The  equations  derived  above  are  then  applied  to 

this  plasma  region  as  a  whole  with  the  radial  dependence  of  the  plasma 

variables  either  assumed  a  priori  or  integrated  out.  For  instance,  one  could 

begin  by  specifying  a  density  and  velocity  profile  consistent  with  the  mass 

continuity  equation.  These  profiles  would  depend  on  a  Lagrangian  coordinate 

spanning  the  plasma  shell,  the  inner  (R.  )  and  outer  (R  shell  radii,  and 
"  in  '  out 

the  inner  (V.  )  and  outer  (V  .)  shell  velocities.  The  momentum  equation  is 
in  out  ^ 

next  integrated  over  the  shell  volume  to  obtain  an  equation  for  V 

Likewise,  the  internal  energy  equation  is  integrated  over  the  shell  volume  to 

obtain  an  expression  for  the  average  internal  energy.  This  approach  was 

2  5 

employed  by  Guillory  and  Terry  ,  McDonald  and  Ottinger  ,  however,  found  that 

the  total  plasma  plus  field  energy  was  not  conserv/ed  in  the  Guillory  and  Terry 

model.  Instead,  McDonald  and  Ottinger  employed  the  kinetic  energy  equation  to 

determine  V  with  the  assumption  V.  =  V  Further,  within  the  plasma 

shell  they  considered  two  subzones  which  differed  in  temperature  but  not 

density:  an  outer  zone  from  R  to  R  ^  of  fixed  mass  which  contained  a 

out  out 

uniform  current  density,  and  an  inner  zone  form  R.  to  (R  ^-S)  which  was 

in'  out  ' 

current  free. 

An  alternative  approach  is  to  view  a  zero-D  model  as  a  special  case  of  a 
typical  full  1-D  code  with  many  zones.  In  this  approach  no  plasma  profiles 
are  assumed  a  priori,  instead  (i)  the  fluid  is  divided  into  a  number  of  zones 
such  that  each  zone  boundary  moves  with  the  fluid  velocity  at  that  interface, 
and  (ii)  spatially  averaged  values  for  the  plasma  variables  are  used  in  each 
zone.  Ihe  definition  of  the  averaged  quantities  are  determined  by 
approximations  to  the  volume  integrals  of  the  fluid  equations.  Since  the 
integration  is  over  a  moving  zone  it  is  appropriate  to  adopt  a  Lagrangian 
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framework  for  the  fluid  equations.  This  requires  the  transformation  of  the 

independent  variable  from  the  spatial  coordinate  r  to  an  indexing  coordinate 

j,  viierein  each  fluid  element  moves  about  between  fixed  values  of  j.  The 

spatial  position  of  a  fluid  interface  with  index  j  is  R(j)=R.,  and  its 

elocity  ts  DR./Dt  -  V..  The  Z-axis  is  at  1=0  and  the  outermost  interface  is 
J  J 

at  j=NJ.  The  Jacobian  for  the  transformation  between  the  two  coordinate  is  9 
3R/9j,  giving  for  the  volume  element  about  j,  2ri/iZRdR  =  2itAZR(  j  )0(  j  )dj  .  To 
integrate  the  dynamic  equations  we  regularly  employ  the  theorem, 


R(j+1) 

Iy  2nR0dj 

R(i) 


Dt 


R(jrl) 

I  fZnrdr 

R(j) 


-|R(jrl) 

^R(j) 


(14) 


where  f  is  an  arbitrary  quantity.  We  note  that  the  meaning  of  Df/Dt  is  that 
the  function  f  is  first  to  be  written  as  a  function  of  the  Lagrangian 
coordinate  j  and  the  time  t  before  differentiating.  The  relation  above  is 
then  readily  proven  by  transforming  the  integral  on  the  right  to  a  function  of 
j  (2nrdr  -*  2riP0dj),  bringing  in  the  time  derivative  and  finally  using  the 
standard  relation  between  the  Lagrangian  and  Eulerian  time  derivative: 


2i  _  li 

Dt  ”  at  0  aj  • 

We  now  integrate  the  continuity,  momentum,  and  energy  equations  over  a 
volume  element  2riR9dj  and  apply  eqn.(14).  For  the  continuity  equation  (1)  one 
finds 


R(j+1)  R(j+1) 

^  j  p.2rtrdr  +  [2TiRp,(Vj^-  §|-)]  =  0,  (15) 

R(j)  R(j) 

which,  is  the  simple  fact  that  the  total  mass  between  Lagrangian  coordinates  j 
and  j+l  is  constant  since  the  grid  velocity,  DR/Dt ,  equals  the  mass  velocity, 
'/j^.  In  tile  volume  integral  of  the  total  momentum  equation  (2)  we  use 
different  upper  a:.d  lower  bcunds: 


8 


R(j+l/2) 


j  +  1/2 


57  J  *  I  ■ 

R(J-l/2)  j-1/2 

For  the  ion  internal  energy  equation  (3)  we  use 
continuity  equation: 


R(j+l/2) 

-  J  ij^B^Ziirdr.  (16) 

R(j-l/2) 

the  same  limits  as  in  the 


R(j+1)  j  +  1  R(j  +  1) 

„,Ci2.rdr»J  (p.*Q„.^)|j(RVg)2ndj  *  [2nRqj  „]  .  J  Qj^2nrdr.  (17) 

R(j)  j  R(j) 

Likewise,  the  volume  integral  of  the  electron  internal  energy  (6)  is 

R(j+1)  j+1  R(j+1) 

57  J  »^t^2nrdr  .j  |j(RVj,)2Tldj  *  [2nRq^  J 
R<J)  i 

R(iti) 

=  J  [a  -  Qig  +  Vz  .  (18) 

R(j) 

Again  for  later  reference  we  list  the  volume  integral  of  the  kinetic 
energy  equation  (11), 


D_ 

Dt 


R(j+l/2) 

J  K''r^^  * 

R(j-l/2) 


(Pi»  Pe» 


R(j+l/2) 

R(j-l/2) 


J 


R(j  +  l/2) 

<Pi"  Pe"  Qvis>?  2nrdr  = 


R(j-l/2) 


R(j+l/2) 

R(j-l/2) 


the  integral  of  the  plasma  energy  equation  (12)  between  arbitrary  boundaries 
R(j=a)  and  R(j=b), 


R(j=b) 

R(j=a) 


P  + 


VIS 


R(j=b) 

R(j=a) 


D_ 

Dt 


r  -,R(j=b) 

A'  V  J  “  /  y 

r  /' 

and  the  integral  of  the  field  energy  equation  (13)  with  the  same  boundaries, 


R(j=b)  g2 


If  I  -  [(tI) 

R(  j=a) 


b:.  ^  -,R(j=b) 

2iiRV_+  ^„B.2jiR 

“  ^  ♦  4(J.a) 


R(J-b) 


J 


J2E22jirdr.  (21) 


R(j=a) 


C.  Finite  Difference  Form  of  the  Equations. 


The  basic  set  of  equations  which  are  used  to  advanced  the  density, 
velocity,  and  internal  energies  are  eqns. (15)-(18) .  We  next  convert  this  set 
from  time  differential  and  spatially  integrated  forms  to  true  finite 
difference  forms.  The  mass  continuity  equation  (15)  determines  the  advanced 
average  density  of  each  zone: 


(n+1) 

Pi,j+l/2 


i+1/2 


(t=0) 


(2,  ) 


where  M.  .  is  the  total  mass  initially  between 
j  +  l/Z 

time  at  step  n. 

For  the  total  momentum  equation  we  take 
velocity  at  interface  R.  between  R.  ,  and  R.  , 

1  j+1/2  j-1/2 

as 


R.  and  R.  -,  and  t*^  is  the 
1  J+1 

to  be  the  mass  weighted 
and  approximate  eqn.  (16) 


+  M  y(r>  +  l) 

.i^l/2^  1-1/2)  ^ 


V) 


(n) 


2AZ 


■J 


At 


[“’j. 


1/2  ^j-1/2^^'® 


j- 


(n+1/2) 


where 


R(j+l/2) 

=  '  [  I 

R(j-l/2) 


(n+1/2) 

,  (23) 
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P  =  p.  +  p  +  Q  . 

In  this  approximation  we  have  associated  half  of  the  total  mass  in  each  of  the 
adjacent  zones  with  the  velocity  at  the  interface  j.  '’he  integral  over  the 
pressure  force  is  evaluated  for  a  step  function  approxiriation  to  the  pressure 
profile  across  the  interface.  The  temporal  finite  differencing  is  indicated 
by  the  superscripts  (n+1)  and  (n).  The  superscript  (n+1/2)  does  not  denote 
the  time  average  of  a  quantity,  but  rather  a  value  determined  from  a  half-step 
predictor . 

The  form  of  the  artificial  viscosity  is 


)  .  .  1 -o  =  b  .  p.  .  1 /oC  .  1,0  min{0,V.-  V. 

vis,j+l/2  VIS  i,j+l/2  s,j+l/2  j  j+1^ 


where  b  .  is  a  constant  parameter  chosen  to  match  the  results  from  a  standard 

VIS 

multi-zone  hydrodynamic  code.  is  the  isothermal  sound  speed  of  the  ions. 
Note  that  the  artificial  viscous  pressure  only  contributes  if  the  zone  is 
undergoing  compression. 

For  eqns.(17)  and  (18)  we  take  the  ion  and  electron  internal  energies  as 
uniform  between  and  Equation  (17)  is  approximated  as 


M.  , 
■1+1/2 

AZ 


(n+1)  _  (n) 

i.i+1/2  1, 1+1/2 

At 


+ 


r'®j+i^i,j+r 


2nR.q.  .1 
1^1. jJ 


(n+1) 


(n+1/2) 

j+1/2 


(24) 


where  is  to  be  evaluated  using  the  spatially  averaged  temperature  and 

density  according  to  eqn.(5).  The  pressure  integral  of  eqn.(17)  has  obviously 

been  approximated  by  extracting  the  averaged  pressure  of  zone  j+1/2  and  then 

carrying  out  the  integral.  The  temporal  location  of  the  heat  flux  is  at  the 

advanced  time,  for  in  practice  one  time  splits  off  this  term  and  calculates 

its  effect  implicitly  after  advancing  the  energy  equation  to  t*^^  .  We  have 

also  not  yet  specified  the  time  location  of  V'.  and  V'  .  in  the  second  term  of 

J  j+1 

this  equation,  but  will  do  so  shortly  by  an  energy  conservation  relation.  The 
corresponding  result  for  the  electron  internal  energy  equation  (18)  is 


11 


(n+1)  _  (n) 

•  i/o  • 


^  (n.1/2)  ^^^.1/2)^,  _  2nR(n-l/2)vO 

AZ  At  *^e,j  +  l/2  I  j  +  1  j  +  1  j  jJ 


[21TR.  ,q  .i-2rtR.q  .1 
L  j+l^e,j+l  J  e,jJ 


(n+1) 


r  _  _  voii r  r 


R(j+1) 

2 

qj^J22iirdr 

R(j) 


T(n+l/2) 

J  ,  (25) 


where  ^(*^■^^■^2)  average  heating  (cooling  if  negative)  rate  over  the 

timestep.  The  same  statements  concerning  the  time  location  of  V'.  V'.  ,  ,  and 

J  J+1 

the  heat  flux  terms  made  for  the  ion  internal  energy  equation  apply  to  this 
equation.  Note  that  the  total  number  of  electrons  in  zone  j+1/2  has  been 
assumed  to  be  a  constant  in  this  approximation.  This  reflects  the  time 
splitting  between  ionization  and  hydrodynamics  wherein  the  number  of  free 
electrons  is  unchanged  during  a  hydro-step.  This  approach  requires  that  the 
total  electron  internal  energy,  (3/2)M  k_T  /m  +  Me  ,  remain  unchanged  during 

6  D  6  6 

an  ionization  step. 

The  approximation  to  the  integrals  made  in  the  above  equations  are  not 
arbitrary,  but  are  determined  by  plasma  energy  conservation  over  the  entire 
system  within  the  finite  difference  scheme.  We  now  derive  through  algebraic 
manipulations  a  plasma  energy  conservation  relation  starting  from  the  finite 
difference  forms  used  to  advanced  the  velocity  [eqn.(23)]  and  the  internal 
energies  [eqns.(24)  and  (25)].  The  result  will  then  be  compared  with  the 
exact  conservation  equation  (20)  as  a  check,  on  the  finite  difference  forms 
(23)-(25).  Define 


<V.> 

J 


v/"))  . 


(26) 


and  likewise  for  <V .  i,~>.  This 

J+1/2  .  T/nx 

•1  ..u  ,7  (n+1/2) 

necessarily  the  same  as  Vj  , 

predictor.  Next  multiply  the  finite 

equation  (23)  by  to  obtain 


time  averaged  velocity 
which  is  the  result 
difference  form  of  the 


at  R.  is  not 
J 

of  a  half-step 
total  momentum 
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1 

2 


r- 


i^l/2 


+  M. 


2AZ 


1/2)  1  r 

J  At  L 

rj  j  ' 

r  j  j  J 

f  (n+1/2) 

n+1/2 


<V.> 


R(j+l/2) 

[J 

R(j-l/2) 


■J2B^2iirdr  I 


,(n-Hl/2) 


(28) 


This  relation  is  the  derived  finite  difference  form  of  the  kinetic  energy 
equation.  If  we  note  that 


/  (n+1/2)  p(n+l/2)) 

1  j+1/2  j-1/2  j 

'  J  J 

I-  ■'R(j-l/2) 

p(n+l/2) 

"  j+1/2 

f2nR("+l/2) 

r  j+1/2  j+1/2 

p(n+l/2) 

■  j-1/2 

f2nR^"'^^^^^<V.>  - 
V  J  J 

J-1/2  ^''j-l/2^j  ’ 

then  the  form  of  the 

terms  on  the 

left  hand  side 

of  equation  (28)  is 

consistent  with  those 

same  terms  in 

equation  (19). 

Equation  (28)  thus 

expresses  the  local  conservation  of  kinetic  energy  over  the  domain  Rj_i/2 
^j  +  l/2‘  other  hand,  the  sum  of  eqns.  (24)  and  (25)  expresses  the  local 

conservation  of  total  internal  energy  over  the  domain  R^  to  Rj+j*  Since  the 
two  domains  are  different  we  cannot  derive  a  finite  difference  form  for  local 
conservation  of  total  plasma  energy  (eqn.  20)  starting  from  the  finite 
differenced  momentum  and  internal  energy  equations.  However,  total  plasma 
energy  can  be  shown  to  be  globally  conserved.  To  see  this  sum  equation  (28) 
over  the  interfaces  between  the  axis  and  the  outer  boundary  (0  <  j  <  NJ)  and 
rearrange  the  result  to  read 


At 


NJ-1 

■  I 

j=0 

NJ-1 

I 

j=0 


1  ^i^l/2 

2  AZ 


p(n+l/2)f 

J+1/2  [ 


1 


2 


2  HR 


^n+l/2)<v.  >  -  2llR^"'^^^^^<V.>  1 

J  +  1  J  +  1  J  J  ) 


2 


13 


NJ-1 

=  I 

j=0 


R(j+l/2) 


<V.> 

J 


[J 

R(j-l/2) 


-J„B .2nrdr 
c  Z  9 


i(n+l/2) 


(29) 


as  long  as  V(j=0)  =  V(j=NJ)  =  0.  Next  sum  the  internal  energy  equations  (24) 

and  (25)  from  j=0  to  j=NJ-l  and  add  to  equation  (29).  We  now  identify  in 

eqns.(24)  and  (25)  with  <V.>.  At  the  axis  one  must  have  q.  „  =  q  =0,  and 

J  1 ,  n  e ,  R 

we  will  assume  that  there  is  no  energy  loss  due  to  conduction  through  the 
outer  wall  at  j  =  NJ.  Then  the  result  is 


At 


NJ-1 

■  I 

j=0 


2 


2 


+ 


N  1^1  M 

^^1/2  L(n^l) 

At  AZ  ri,j  +  l/2 


e 


(n) 

i ,  j  +  1/2 


+ 


.(n+1) 
e, j+1/2 


(n)  ) 

e, j+1/2  J 


+ 

+  1/2 


R(j+1) 

J  '1j^'l2^2Tirdr 

R(j) 


R(j+l/2) 

<V.>r  -J„B.2jirdr 
J  J  c  Z  <p 

R(j-l/2) 


(n+1/2) 

.(30) 


This  expression  is  the  derived  finite  difference  form  of  the  total  plasma 
energy  equation  for  the  entire  domain.  Note  the  similarity  of  this  result 
with  the  sum  over  all  the  zones  of  the  exact  expression,  eqn.(20). 

Equation  (30)  is  actually  quite  remarkable  for  the  case  of  no  magnetic 
fields.  It  states  that  one  can  transport  internal  energies  in  a  Lagrangian 
code  and  yet  maintain  global  conservation  for  the  total  energy  (internal  + 
kinetic).  This  result  is  not  based  upon  an  expansion  for  a  small  timestep, 
but  rather  arises  from  the  finite  difference  forms  (23)-(25)  for  the  momentum 
and  internal  energy  equations.  The  proper  sequence  of  steps  to  achieve  energy 
conservation  in  a  hydrodynamic  system  is  as  follows; 


(s-i)  Using  the  values  at  t^"^^  do  a  half-step  predictor  to  find  p. ,  V, 

.  .n+1/2  f  .  .  ^ 

E.  ,  £  at  t  for  each  j . 

(s-ii)  From  these  quantities  compute  p^,  p^,  and  R  at  at  each 

j.  (If  there  is  a  magnetic  field,  also  compute  the  currents  I  and 
J„  and  the  field  B.  at  t*^^  as  discussed  below.) 

(s-iii)  Compute  '  at  each  j  from  eqn.(23). 
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(s-iv)  Form  from  eqn.(26). 

(s-v)  Use  these  quantities  in  the  energy  eqns.(24)  and  (25)  to  find  the 
final  advanced  quantities  and  at  each  j. 

Unfortunately,  with  the  addition  of  magnetic  fields  and  an  external 
circuit,  conservation  of  the  total  plasma  plus  field  energy  is  not  as  straight 
forward  as  it  is  for  a  hydrodynamic  system.  We  have  evaluated  the  J>a  force 
term  in  eqn.(23)  and  the  resistive  heating  term  in  eqn.(25)  at  but 

have  yet  to  show  consistency  with  the  derived  plasma  plus  field  energy 
conservation.  We  will  address  this  problem  in  §11. G. 


D.  Hydrodynamics  for  the  Zero-D  Model. 

One  of  the  objectives  in  deriving  the  above  general  equations  is  to 
explicitly  demonstrate  the  underlying  assumptions  and  limitations  of  a  zero-D 
dynamic  model.  To  reduce  the  above  set  of  equations  to  two  momentum  and  two 
internal  energy  equations,  begin  by  considering  three  zones  for  the  annular 
shell  plasma.  The  inner  zone  extends  from  the  axis  at  j=0  to  j=l  (R^),  the 
central  zone  from  j=l  (R^)  to  j=2  (R2)>  and  the  outer  zone  from  j=2  (R2)  to 

j=3  (R^).  This  last  radius  is  that  of  a  fixed  cylindrical  boundary  termed  the 
wall.  The  following  assumptions  are  made; 

(a-i)  The  total  mass  of  the  central  zone  (j  =  1+1/2)  dominates  the  mass 
of  the  inner  and  outer  zones,  ^-^^2  ^3/2  ^5/2’ 

(a-ii)  The  temperatures  of  all  three  zones  are  nearly  equal. 

(a-iii)  The  ion  and  electron  number  densities  of  the  central  zone  are 
always  much  larger  than  those  of  the  inner  and  outer  zones. 

(a-iv)  The  outer  zone  is  current  free:  =  0  and  =  2I/rc  for  r  >  R2. 

Since  (a-ii)  and  (a-iii)  also  imply  that  the  viscous 

pressure  of  the  central  zone  likewise  dominates.  For  an  ideal  gas  law  the  same 
condition  holds  for  the  thermal  pressures  and  the  total  internal  energies  (Me) 
of  the  central  zone.  Under  these  assumptions  the  momentum  equations  from 
eqn.(23)  for  and  V2  can  be  written  as 
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3/2  ^  U  ^  (n+l/2)„  (n+1/2) 

2AZ  At  3/2  1 


R(j=3/2) 

[  J  ^Jz®^2»irdr] 

R(j=l/2) 


(n+1/2) 


(31a) 


M  R(j=2)  fn+1/2) 

"3/2  2  ''2  (n+l/2)„  _(n+l/2)  ff  lTn'’nM.r  ^  mim 

2AZ  Ft  -  ^3/2  ^’^2  =  ~  U 

R(j=3/2) 

The  upper  limit  of  this  last  integral  stops  at  R2  according  to  (a-iv).  To 
eliminate  the  heat  flux  terms  in  the  internal  energy  equations  we  separately 
sum  eqns.(24)  and  (25)  over  all  three  zones.  Assumptions  (i)  and  (ii)  allow 
us  to  write  these  sums  as 


„  g(n+l)  _  (n) 

3/2  1,3/2  1,3/2 


fr.  n  \("+l'^2)  f-  p(n+l/2)  .  -  _(n+l/2) 

"  (Pi"  Qvis>3/2  n  <^^2^  - 


[g  Vol](""^^^^  .32_v 

=  [*^16  ,,,  ’  <^2a) 


j=3/2 


M  -  £("^ 

e,3/2  e,3/2  e,3/2 


+  p("3}2^^  [2iir(""^^^^<V2>  -  2rtR(""^^^^<V^>j 


=  [<^  -  =u>  '^] 


(n+1/2) 

j=3/2 


[J  -'I 


J2^2itrdr  j 


(n+1/2) 


(32b) 


The  radiative  term  A  only  contributes  from  the  central  zone  since  it  has  the 

dominant  number  of  radiators.  The  ion-electron  energy  exchange  term  only 

2 

contributes  from  the  central  zone  since  Q.  ®  n  . 

le  e 

The  assumption  (a-ii)  is  reasonable  as  long  :.s  the  heat  transfer, 

particularly  due  to  the  ion  heat  flux,  is  large  compared  to  the  compressional 

2 

heating:  3.9t^^ A(kgT^/mj)/Ar  >>  AV/Ar.  However,  assumptions  (a-ii)  and  (a- 
iii)  do  break  down  when  the  plasma  shell  bounces  off  the  axis.  The  pressure 
of  the  inner  zone  must  then  rise  above  that  of  the  central  zone  as  the 
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collapse  is  reversed.  Further,  if  there  is  some  anomolous  resistivity  in  the 
outer  lov  density  zone,  then  assumption  (a-iv)  may  also  be  violated. 

Equations  (22),  and  (31a)-(32b)  are  the  reduced  set  of  equations  for 
advancing  the  hydrodynamic  quantities  of  an  imploding  shell  in  our  zero-D 
code.  In  the  case  of  a  filled  plasma  cylinder  the  interface  at  j=l  is 
identified  with  the  Z-axis.  The  momentum  equation  (31a)  for  Vj  is  no  longer 
needed  and  the  internal  energy  equations  determine  and  in  the  region 
between  the  Z-axis  and  R2.  The  problem  with  assumptions  (a-ii)  and  (a-iii)  at 
the  time  of  the  bounce  has  essentially  been  integrated  away. 


E.  Magnetic  Induction  Equation. 


One  of  the  most  important  dynamical  effects  in  a  Z-pinch  is  the  amount 
of  resistive  heating  inside  the  plasma  during  the  run-in  phase.  This  is  a 
consequence  of  the  diffusion  of  the  magnetic  field  from  the  outer  plasma 
boundary  into  the  collapsing  plasma  shell.  Rather  than  fix  a  skin  depth  or 
define  a  current  density  profile  in  the  plasma  shell,  we  will  develop  an 
approximate  solution  for  the  magnetic  diffusion  in  the  moving  plasma.  The 
induction  equation  describing  the  diffusion  is  found  by  substituting  the 
expression  for  E^  from  Ohm's  law,  eqn.(8),  into  Faraday's  law,  eqn.(9),  and 
using  the  expression  for  from  eqn.  (10): 


3B  . 
—1 
at 


(33) 


where  we  have  dropped  the  thermal  friction  term  in  for  the  reason  stated 
below  eqn. (13).  To  solve  eqn. (33)  we  make  several  additional  assumptions  to 
those  listed  in  §11. D; 


(a-v)  The  current  skin  depth  (5)  inside  R2  is  sufficiently  small  that 

B,/r  «  aB./ar. 

9  9 

(a-vi)  The  resistivity  is  spatially  uniform  over  the  diffusing  region, 

=  n^(t). 
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To  convert  this  equation  from  an  Eulerian  framework  to  a  Lagrangian  one 
note  that 

a  3®*  3B  .  B  3V_  DB .  B .  ^  DeB . 

_ i.  .  i-(v  B)  -  — S+V  _  1  t 

3t  3r^  "  3t  R  3r  e  3j  "  Dt  9  Dt  "  0  Dt 

Finally,  note  the  specific  relation 

R(j,t)  =  R^(t)  +  (R2(t)  -  R^(t))j  =»  e  =  ||  =  R2(t)  -  Rj(t)  =  M.  (34) 

Equation  (33)  can  now  be  written  as 


D0B 

“dF 


2  2 
c  n,  3  0B. 

_ i  P 

2  2 

4 110  33 


(35) 


Note  that  if  =  0  for  t  >  tj,  this  equation  states  that  the  magnetic  flux 
AZARB^  remains  fixed  at  its  t^  value  even  though  6R  may  continue  to  change. 
Next  transform  to  a  new  time  variable; 


■  I 


t  c^nj^(t') 


dt ' 


'0  4ll0  (t') 

and  thereby  reduce  eqn.(33)  to  a  planar  diffusion  equation: 


D(0BJ  3^(eBJ 

_  t  =  +  _  .T 


Dt 


9j 


.2 


The  solution  is  subject  to  the  initial  and  boundary  conditions 


(36) 


(37) 


0B .  =  0  at  T  =0  and 
9 


e(T) 


21  (j) 

R2('i:)c 


at  j  =  2  (r=R2), 


(38) 


The  boundary  condition  at  the  outer  plasma  radius  reflects  the  value  of  the 
magnetic  field  in  the  outer  vacuum  adjacent  to  the  plasma  due  to  the  current  I 
within  the  plasma.  Equation  (37)  can  be  solved  by  the  method  of  Laplace 
transforms.  The  result  for  the  magnetic  flux  is 
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(2-i)  /4(t  t' ) 


For  the  following  manipulations  set 


♦(  T)  =  0(  T) 


(40) 


R2(t)c  • 

In  order  to  make  the  solution  (39)  tractable  for  numerical  computation  we  will 

;  i  to 
n+1/2 


take  #  to  be  constant  over  each  time  interval  ( t''*,  t"^^^) ,  and  denote  its  value 


at  T  by  ^rj^j/2'  Equation  (39)  can  then  be  written  as 


0(t)B^(t) 


N+1/2 


J| 


t  ^-(2-j) 


dx' 


f') 

*  n^O  ^  *"*''2 


I, 


x"*^  2 

.-(2-j)V4(x-x') 


dx' 


Next  making  the  change  of  variables  from  x'  to  =  (2-j )/[ 2/(x-x' ) ] ,  the 
solution  can  be  calculated  at  time  >  x  >  x*^  from 


where  erf(x)  is  the  error  function  and  s  2/(x-x*^).  Note  that  Oq  =  2/x  > 

>  •••  > 

One  can  interpret  this  last  equation  by  noting  that  at  the  time  step 
^n+1/2’  ^  component  of  the  magnetic  flux  ^^+1/2  plasma  and  the 

diffusion  length  scale  from  that  time  on  is  R2“n’  at  the  outer  radius 

where  j  =  2  the  boundary  condition  eqn.(38)  is  satisfied,  and  at  the  initial 
instant  when  aQ  =  0  there  is  no  magnetic  flux  in  the  plasma. 

For  the  current  density  AnJ^/c  =  B^/R  +  (3B^/9j)/0  one  can  calculate 
directly 
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c 

4tx 


O  3/2„2 
2  ri  e 


r  N 


-(2-j)^/o  ^ 

V  J  /  p, 


n+1/2  '^n-1/2^ 


2  2 
-(2-j)Va 

»  *1/2  - 5 - ^  ]•  <‘2) 

o 

A  maximum  for  the  skin  depth  6  of  the  current  in  the  plasma  shell  at  "time"  x 
can  be  determined  by  the  penetration  depth  of  the  first  component  of  the 

magnetic  flux  *^^2  enter  the  plasma  at  time  T=0:  S  ~  a^fR^-R^^).  This  value 
is  clearly  an  overestimate  of  the  skin  depth  since  the  magnetic  flux 
components  which  enters  the  plasma  shell  at  later  times,  *n+i/2’ 

significantly  larger  and  have  a  smaller  penetration  depth. 

The  total  current  within  the  radius  R~  at  a  time  x  lying  in  the  interval 

ki  ki  i 


I(T)  = 


J2(2-j,x)2fl(R2(x)-(2-j)0(x)]0dj 


In  the  limit  of  a  narrow  skin  depth  one  can  neglect  the  second  term  in  the 
square  brackets  in  the  integrand  and  the  lower  limit  of  the  integral  can  be 
replaced  by  One  then  finds  by  carrying  out  the  integration  the  consistent 


‘■N+1/2'  Under  the  same  estimates  for  computing  the  integral, 
the  resistive  heating  at  time  in  the  electron  internal  energy  equation 


result  I(x)  =  Ij 
the  resistive  h( 
(32b)  is  explicitly 


n2 


2nR0dj 


c^R2 


1/2 


^*1/2^ 


[  T 


N  N 

^  I  Z  (*, 


^/*n+l/2  *n-l/2^  a 
n=l  n 


2/2 


,  ,  n  +  1/2  '^n-l/2^^'^m+l/2  ^m-1/2^  .  „  „;i/2 

n=l  m=l  f  2  2 

n  ^m  a  4  a 

(  n  m 


n^l 


l/2'"n+l/2  *n-l/2^ 


2/2 


(  2  2 
a  +  a 
(  n  o 


1/2 


(43) 
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Likewise  the  mean  value  for  the  magnetic  field  within  is 

^  f  K^1/2“n  ^  Jq  Vl/2(“n+l-“n>]  ' 

The  Maxwell  stresses  of  eqn.  (31a)  and  (31b)  can  be  readily  evaluated  in  this 
limi t : 


2 

=  -  -  B^<j=3/2)]  =  -  ^  ^  (45) 


and 


c  Z  <1 


)2itR0dj 


-  K^J(:i=3/2)  . 


(46) 


F.  Circuit  Model 


The  electromagnetic  force  driving  the  current  which  implodes  the 

plasma  shell  is  provided  by  a  generator  external  to  the  plasma  load.  The 

entire  system,  plasma  plus  generator  or  machine,  can  be  modelled  by  a  simple 

circuit.  Let  the  machine  have  a  constant  resistance  Z  and  inductance  L  , 

m  m 

while  the  plasma  sustains  a  varying  voltage  drop  at  the  outer  boundary  or 
wall  (r  =  R^).  The  equation  for  this  circuit  is 


L  ^ 
m  dt 


Z  I  = 
m 


V  (t) 
w'  ' 


(47) 


To  find  V^,  integrate  Faraday's  law,  eqn. (9),  over  the  region  between 
the  outer  plasma  radius  R2  and  the  wall  radius  R^.  From  relation  (14) 
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R(j=3) 

^  I  B*'**  -  *  <^*>3.2  -  '  <^^,(82)  . 

R(j=2) 

Since  Vn(j=3)  =  0  we  have 
K 

R(j=3)  y 

=  E2(R3)AZ  =  ie^dRAZ  .  AZ(B^-|  ).^2  "  ‘ 

R(j=2) 

According  to  the  discussion  following  eqn.(13),  we  will  ignore  the  term  in 
the  expression  for  from  Ohm's  law,  eqn.(8).  The  last  two  terms  of  the 
above  equation  then  reduce  to  AZhj^C  t  )J2(  j=2 ,  t )  .  Since  =  2I/Rc  in  the  outer 
region  [assumption  (a-iv)  of  §11.0]  the  integral  in  eqn.(48)  can  be  done 
explicitly  and  the  time  differentiation  carried  out.  The  result  for  the 
potential  drop  at  the  wall  is 


where 


V 

w 


— 
'geo  dt 


Z  I 


geo 


and 


Z 

geo 


(49) 


(50) 


are  the  geometrical  inductance  and  motional  resistance  arising  from  the  plasma 
load.  Using  eqn.(50)  in  eqn.(47)  gives  the  working  equation 


dl 


(L  +  L  )  ^  +  (Z  +  Z  )I 

in  geo  dt  m  geo 


=  ^mf  “  n^J2(R22t)AZ.  (51) 

This  equation  is  solved  using  a  finite  difference  scheme: 


[ (n.1/2) 
geo 


j(n>l)_  j(n) 

Ft 


.(n+l/Z) 

geo 


l(n+l)^  l(n) 
2 


n  +  1/2 


>mf  -  Vz(''2>^]  •  (^2) 
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F.  Total  Energy  Conservation. 


The  plasma  plus  field  total  energy  equation  is  the  sum  of  eqns.(12)  and 
(13).  The  integral  of  the  result  from  the  axis  to  the  outer  boundary  is  the 
sum  of  eqns.(20)  and  (21).  If  there  is  no  heat  flux  through  the  outer  wall 
and  eqn.(8)  is  used  for  this  later  sum  is 


D_ 

Dt 


R(j=3) 

I  (“I'i* 


0 


+ 


D_ 

Dt 


2  nr  dr 


0 


R(j=3) 

=  +  J  A2nrdr,  (53) 
0 

where  is  the  potential  drop  at  the  wall  according  the  eqn.(48).  Note  that 
IV^  is  the  input  power  to  the  plasma  load  from  the  machine  while  the  second 
term  on  the  right  hand  side  of  eqn.(53)  is  the  negative  of  the  radiative  loss 
rate  per  unit  length. 

The  objective  of  this  section  is  to  derive  a  corresponding  expression 
for  eqn,(53)  starting  from  the  difference  equations  used  to  advance  V,  e^,  e^, 
and  V^,  and  the  differential  equation  for  B^.  Part  of  this  procedure  has 
already  been  done;  the  total  plasma  energy  equation  was  derived  in  eqn.(30). 
To  derive  an  energy  equation  for  the  field  let  us  start  with  the  expression 
for  the  energy  input  from  the  machine  over  the  time  At.  From  eqn.(49) 


At<I>V 


(n+1/2) 


w 


=  At 


(n) 


)K 


n+l/2) 

geo 


j(n+l)_  j(n) 
At 


(n+1/2) 

'geo  2 


f  x(n+l/2), 

+  Az(n^J2(R2))  J  .  (54) 


Using  the  definitions  in  eqn.(50),  we  can  rearrange  the  right  hand  side  of 
this  equation, 


,(n+l).2  R  V  At, (n+1/2) 

( V-  J  ^  huj)  - 


(n+1/2) 
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,(n+l)_(n)  f.V„At,(n+l/2) 

^  “  [irJ  * 


P  ,(n+l/2 

[Vz(R2)J  AZAt. 


(55) 


Note  that 


r  /^3'(  -.(n+1/2)  ,  f  ^2  ) 

]  =  ^"(j^(n+l))  ^  (j^(n+l/2)  ) 


(n+1) 


v(n.l/2)^^ 


r  ^3  )  (  )  r  <v2>At>^ 

=  ^"(^(n+lTJ  ^  lj^(n+l/2)  j  +  2<R2>  J  ~  ^^“OTTlTiy  ' 

2  w  2 

As  long  as  |V2At/R2|  <<  1,  the  sum  of  the  last  three  terms  is  much  less  than 

the  first  term,  which  grows  as  the  pinch  implodes.  A  similar  expansion  can  be 

obtained  for  the  second  term  in  expression  (55).  Further,  we  note  that  for  R2 

<  r  <  R-,,  B,  =  2I/rc,  which  implies 
3  9 

R(j=3)  .2  0  p 

I  • 

R(j=2) 

Using  the  above  relations  in  eqn.(54)  gives 


R(j=3)  „2 

[  J  ( 

R(j=2) 


(n+1) 


R(3=3)  j,2 

[  J  (  8^)^--^] 

R(j=2) 


(n) 


/  /  \  17(0+1/2).. 

.1  (n+1/2)  i(o+l)j(n)  V2  At 

=  "  - 2 - „(n+l/2")  - 

C  ^<1 


r  , (n+1/2) 

[Vz(R2>J  • 


We  next  work  on  the  magnetic  field  energy  in  the  region  r  <  R2. 


R(j=2)  g2 


[  J  (  8^)2'"^'“^] 
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R(j=2)  2  .  . 

[J  (a-JM 
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(n.l)B(n.l)  ^  _  e(n)B(n))dj 


+ 
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l^(n+l)g(n+l)B(n)g(n)|'  9^ 

A  <|)  <l>  I  ^(n) 


For  small  timesteps 


0(n+l) 


(57)  , 

! 
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Q(n+l)B(n+l) 

♦ 


e(n)B(n) 


=  At 


(n+1/2) 


At 


(n+1/2) 


by  eqn.(35),  and 


Q(n)  0(f'+l) 

^(n)  ~  jj(n+i) 


1_  (  <V>^t) 

3j  r  <R>  } 


Using  these  last  tvo  relations  in  eqn.(57)  and  integrating  by  parts  gives 


R(j=2)  b2 


R(j=2)  gZ 


r  r  n(n+l)  r  r"'" 

[  J  ( -  [  I  ( 

0  0 

1  r/  r  -^(11+1/2) 


,(n+l) 


1  /  N  V  ^  ^  (  an)  <V„> 


j=2 

-  "J  I  [(jzH 
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(n+1)  .  >.(n+l/2) 


25 


"I  ^  I  . »  “i- (58) 


Finally  add  eqns.(56)  and  (58),  and  note  that  6^(82)  =  21/R2C  to  find 


0  0 


R(j=3)  .2 

t*  /•  D  ,  s 


^I>,(n.l/2) 

AZ  w 


r  ^(n+i/2) 

J  2  [KH  ^  JlVzJ  2.dj 


M  —  dK)  M  ^  K)  M  (59) 

0 

Equation  (59)  is  the  derived  form  of  the  field  energy  equation  (21).  We 
note  that,  unlike  the  derived  total  energy  equation  for  the  plasma  (eqn.(30)], 
eqn.(59)  is  only  valid  to  first  order  in  At.  The  derived  form  of  the  total 
plasma  plus  field  energy  equation  is  given  by  the  sum  of  eqns.(30)  and  (59). 
Comparison  of  this  result  with  eqn.(53)  indicates  the  proper  time  location  for 
the  integrand  in  the  resistive  heating  and  JxB  work  terms,  the  second  and 
third  terms  on  the  right  hand  side  of  eqn.(59).  However,  the  required  timing 
of  these  two  terms  presents  a  problem.  At  step  (s-iii)  of  §11. C  for  the 
computation  of  the  quantities  and  are 

9  o 

not  known.  But  we  note  that  to  first  order  in  At  we  can  make  the  replacements 


1  U  N(n  +  1).  .(n)  ,  >.(n),  .(n+1)^  .tn+i/^^; 

2-  Lr  J  M  ^  («%]  M  ]  i>  K^z"®) 

Now  the  expressions  for  the  resistive  heating  in  eqns.(30)  and  (59)  agree. 
However,  for  the  JxB  work  term,  <V>  is  inside  the  integral  in  eqn.(59),  but  in 
eqn.(30)  the  integral  is  split  up  over  the  j-domain  and  <V>  is  extracted.  In 


(n+1/2) 


a  full  1-D  code  there  are  many  zones  so  that  the  splitting  of  the  integral  is 
not  a  severe  approximation.  In  a  zero-D  code  the  Maxwell  stress  term  in  the 
momentum  equation  (23)  remains  a  potential  problem  for  the  conservation  of 
total  plasma  plus  field  energy.  To  ameliorate  this  potential  difficulty  we 
have  only  considered  applications  of  the  zero-D  code  to  rapidly  collapsing 
plasmas  so  that  the  subsequent  rapid  resistive  heating  maintains  a  narrow  skin 
depth  for  the  current  density  (see  assumption  (a-v)  of  §II.EJ.  Thus  the 
Maxwell  stresses  contribute  only  over  a  narrow  region  at  the  outer  boundary  of 
the  plasma  where  V  =  ^2'  The  estimates  in  eqns. (43)-(46)  are  to  be  evaluated 
as  part  of  the  half-step  predictor  sequence  in  step  (s-ii)  of  §11. C. 

To  maintain  complete  energy  conservation  jV^At/Rj |  must  be  <<  1  as 
required  in  the  above  derivation  for  field  energy  conservation.  In  practice 
we  have  found  that  total  energy  is  conserved  to  ~2^  over  an  entire  simulation 
with  the  zero-D  code. 


III.  ATOMIC  MODEL  AND  RADIATION  TRANSPORT 

o 

The  Al  atomic  model  is  described  in  detail  by  Duston  and  Davis  .  There 
is  excited  level  structure  in  the  eleven  highest  ionization  stages  (H-like 
through  Na-like);  only  ground  states  are  carried  for  the  other  stages.  For 
each  atomic  level  included  in  the  model,  a  rate  equation  of  the  form 


df . 
_ 1 

dt 


Z  W. .  f .  -  Z  W. .  f . 
J  Ji  1  J  ij  J 


(60) 


is  constructed,  where  fj^  is  the  fractional  population  of  level  i,  and  is 
the  net  reaction  rate  describing  the  transition  from  initial  state  j  to  final 
state  i. 


For  the  plasma  densities  modeled  in  these  studies,  the  effective 
populating  and  depopulating  rates  are  generally  fast  compared  with  the 
hydrodynamic  response.  An  equilibrium  assumption  can  be  made  in  which  the 
explicit  time  dependence  is  dropped;  in  this  case,  the  plasma  is  said  to  be  in 
collisional-radiative  equilibrium  (CRE),  wherein  the  plasma  ionization  state 
responds  instantaneously  to  changes  in  hydrodynamic  parameters.  The  processes 
included  in  the  transition  rates  and  the  methods  used  in  calculating  the 
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corresponding  rate  coefficients  are  described  elsewhere  .  Once  the  rate 
equations,  including  radiation  transport,  have  been  solved,  the  electron 
density  and  the  ionization  and  excitation  energy  in  the  plasma  can  be 
calculated . 

Radiation  emission  from  a  plasma  and  its  opacity  are  dependent  on  the 
local  atomic  level  population  densities.  Except  for  optically  thin  plasmas, 
however,  these  populations  also  depend  on  the  radiation  field  through  optical 
pumping,  photoionization  and  photoexcitation.  Thus,  the  atomic  physics 
processes  and  radiation  transport  are  strongly  coupled  and  must  be  solved 
self-consistently.  In  this  model,  an  iterative  procedure^  is  used  in  which 
the  radiation  field  from  the  previous  iteration  is  used  to  solve  for  new  f^, 
from  which  a  new  radiation  field  is  calculated  and  populations  are 
recalculated  until  convergence  is  achieved. 

A  hybrid  model  is  used  to  transport  the  radiation.  Free-free  radiation 
is  treated  on  a  multifrequency  basis.  The  line  and  bound-free  radiation  is 
transported  via  a  probabilistic  model^^. 


IV.  AVERAGE-ATOM  MODEL 


A  time-dependent  average  atom  model  has  been  developed  which  solves 
hydrogenic  rate  equations  for  a  set  of  orbital  population  levels  P^,  where 
1^<10.  These  rate  equations  are  nonlinear  and  are  of  the  form 


dP 
_ r 

dt 


-P  Z  D(n^)  0 
n  m  m 


0  L  C(m->n)  P 
n  m  m 


-PI  + 
n  n 


Q  R 
n  n 


(61) 


where  Q  =(1-P  /g  ) ,  I  is  the  ionization  rate  of  level  n,  R  is  the 

recombination  rate  into  level  n,  D(n->m)  is  the  destruction  rate  of  level  n  by 

transition  to  level  m,  C(m->n)  is  the  creation  rate  of  level  n  by  transition 

2 

from  level  m,  and  g^  is  the  degeneracy  of  level  n  and  is  given  by  2n  for 

hydrogenic  levels.  The  model  includes  radiative  and  three-body  recombination, 

collisional  excitation  and  de-excitation  and  radiative  decay  in  the  rates. 

12 

This  model  is  based  largely  on  the  work  of  Post,  et.  al. 

Once  the  populations  P^  have  been  calculated,  the  average  charge  and 
radiative  cooling  rate  can  be  determined.  An  approximate  way  of  generating  nl 
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subshell  populations  and  energies  has  been  devised;  in  this  way,  /in=0 
transitions  can  be  treated.  This  treatment  of  An=0  transitions  differs  from 
the  more  complicated  approach  of  Post,  et.  al.  Bremsstrahlung,  bound-free, 
and  line  radiation  are  included  in  the  radiative  cooling. 

The  rates  in  the  Al  atomic  model  used  in  these  studies  were  not 
considered  reliable  for  electron  temperatures  below  about  20  eV,  so  the 
average-atom  model  was  sometimes  used  for  electron  temperature  lower  than 
this . 


V.  APPLICATION  TO  ALUMINUM  VAPORS 
A.  Plasma  Shells  of  Varying  Thickness 

The  ZPIMP  code  described  in  the  above  sections  was  used  to  simulate  the 
implosion  of  aluminum  plasma  shells  of  varying  mass  loading  and  initial  radii. 
Mass  loadings  of  30,  35,  and  40  pgm/cm  were  considered.  Initial  outer  and 
inner  radii  pairs  were  (1.00,  0.90),  (1.05,  0.85),  and  (1.15,  0.75)  cm;  these 
pairings  gave  initial  shell  thickness  of  0.1,  0.2,  and  0.4  cm.  Initial 
electron  and  ion  temperatures  of  5eV  were  assumed.  The  coefficient  of  the 
artificial  viscosity,  was  set  at  0.25  after  comparison  of  some  test  runs 
with  a  multi-zone  1-D  Lagrangian  code.  The  resistivity  is  given  by 
Braginskii ' s”  form  in  the  high  field  limit.  For  the  machine  inductance  we 
used  L^  =  58  nanohenrys;  for  the  machine  impedance  =  2  Ohms;  and  the  open 
circuit  voltage  was  chosen  to  match  GAMBLE  II. 

Figure  1  shows  the  minimum  outer  radius  R,  .  =  min{R„(t)}  obtained  in 
this  series  of  runs  (the  minimum  inner  radius  always  went  to  zero).  The 
smallest  of  the  R„  .  always  occurs  for  the  lowest  mass  loadings  and  thinnest 

^  y  fTi  i  n 

shells.  Figure  2  gives  the  times  of  peak  compression  for  this  set  of 
calculations.  The  earliest  compression  times  again  occur  for  the  lowest  mass 
loadings  and  thinnest  shells. 

This  pattern  of  the  extremum  occurring  for  the  lowest  mass  loadings  and 
thinnest  shells  also  occurs  in  the  next  several  plots  (Figs.  3-7)  for  electron 
temperature,  average  charge  and  radiated  power.  The  higher  radiative  cooling 
under  these  conditions  allows  the  plasma  to  be  compressed  to  higher  densities 
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(or  smaller  radii).  The  exceptions  are  the  line  and  continuum  radiation 
maxima  below  the  K  edge,  which  are  fairly  flat. 

All  the  Z  maxima  are  well  above  11,  which  means  that  the  A1  ions  are 
usually  in  the  He- like  and  H-like  states  at  peak  compression.  This  is 
reflected  in  the  fact  that  the  line  and  continuum  radiation  maxima  above  the  K 
edge  are  higher  than  those  below  the  K  edge.  These  results  indicate  that, 
under  the  set  of  conditions  simulated  by  these  calculations,  it  should  be 
possible  to  create  highly-ionized  Z-pinch  A1  plasmas  in  the  NRL  GAMBLE  II 
device. 

B.  Filled  Cylinder  Plasmas 

Using  the  same  mass  loadings,  a  set  of  filled-shell  (inner  radiussO) 

runs  were  made  for  different  initial  outer  radii  ranging  from  0.5  cm  to  2.0 

cm.  Initial  temperatures  of  5eV  were  assumed.  The  next  several  figures  show 

some  results  from  these  calculations. 

Figure  8  shows  the  minimum  outer  radius  .  obtained  at  each  mass 

Z,min 

loading  vs.  initial  outer  radius  R-  .  Generally,  the  trend  is  for  R„  .  to 

2,0  2, min 

increase  with  R.,  .  However,  the  35  pgm/cm  curve  peaks  slightly  at  R„  =1.875 

Z.  f  O  Z  y  o 

cm,  and  the  40  (jgm/cm  curve  shows  a  very  pronounced  dip  between  1.625  cm  and 
2.0  cm.  Unlike  Figure  1,  the  smaller  R„  .  occurs  at  the  larger  mass 
loadings . 

Figure  9  gives  the  times  of  maximum  implosion  (R„  =  •  )  vs.  R.,  for 

z  z  f  Q 

the  diffeieiiL  mass  loadings.  The  larger  times  occur  at  the  larger  mass 

loading  and  larger  R.,  ,  generally. 

Figure  10  displays  the  maximum  electron  temperatures  achieved  during 

these  calculations.  These  maxima  usually  occurred  at,  or  just  prior  to, 

maximum  implosion.  The  higher  maximum  temperatures  result  at  the  lower 

masses.  The  40  ugm/cm  curve  drops  very  rapidly  in  the  area  where  R2  also 

decreases  rapidly,  tlieu  it  seems  to  stabilize  as  R.,  .  increases  again. 

‘  •  _  2 , min  ® 

The  maximum  average  ionic  charge  Z  is  shown  in  Figure  11.  The  results 

seem  to  peak  aiound  R2  ^=1.0;  the  highest  ionization  occurs  at  the  smallest 

mass.  The  40  ugm/cm  curve  crosses  Z=ll  in  the  region  where  R„  .  drops 
_  ^  2,min 

rapidly.  When  Z  is  11  or  greater,  He-like  and  H-like  ions  are  predominant. 
When  Z  drops  below  11,  L.  shell  radiation  appears.  The  result  of  Z  falling 
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below  11  is  enhanced  L-shell  and  total  radiation  as  will  be  seen  in  subsequent 
figures. 

Figure  12  shows  the  maximum  total  radiated  power  for  these  runs.  The 

higher  radiation  occurs  for  the  larger  mass,  which  helps  to  explain  the  fact 

that  the  smaller  R~  .  occur  with  the  larger  mass. 

2, min  ° 

Figures  13  and  14  give  the  maximum  line  and  continuum  radiation  below 
and  above  the  K  edge.  Figures  12,  13a,  and  14a  shew  an  increase  in  radiated 

power  at  40  ugm/cm  at  R.,  =1.69  cm.  This  is  the  point  at  which  maximum  Z 

^  y  O 

crosses  11  and  at  which  R„  .  decreases  rapidly.  The  enhanced  cooling  of  the 

plasma,  as  shown  in  Figure  10,  reduces  plasma  pressure  and  allows  for 

increased  compression.  Beyond  R,  =1.75  cm,  the  radiative  cooling  drops; 

plasma  thermal  energy  is  retained  and  it  becomes  harder  to  compress  the 

plasma.  Hence  R^  .  increases  again. 

^  2, min  ® 

It  should  also  be  noted  in  Figures  13a  and  14a  that  the  L-shell 
radiative  cooling  reaches  a  minimum  at  large  R„  and  rises  again  for  the  35 

Z  y  o 

ygm/cm  mass.  This  accounts  for  the  peak  in  this  R,  .  curve  noted  earlier. 

z,min 

These  curves  show  high  radiative  losses  at  low  R«  for  all  the  mass 

2  j  o 

loadings;  hence,  all  R^  .  are  small  at  small  R^ 

y.fnin  y.n 


VI.  SUMMARY 


The  equations  and  structure  of  the  new,  zero-D,  ZPIMP  code  have  been 
presented.  This  code  combines  the  best  features  of  previous,  but  separate, 
zero-D  codes,  such  as  the  coupling  of  a  CRE  ionization  model  with  radiation 

3 

transport  in  SIMPLODE  ,  and  a  self-consistent  circuit  equation  as  in  McDonald 
and  Ottinger's^  code.  However,  ZPIMP  uses  a  completely  different  approach  for 
the  hydrodynamics:  the  continuity,  momentum,  and  energy  equations  are  solved 
in  a  finite  difference  form  similar  to  a  standard  multi-zone  Lagrangian  code. 
Furthermore,  rather  than  specifying  a  profile  for  the  current  density,  ZPIMP 
uses  an  analytic  solution  to  the  induction  equation  for  the  diffusion  of  the 
magnetic  field  in  a  moving  medium.  The  form  of  the  total  (plasma  +  field) 
energy  conservation  relation  is  explicitly  derived  from  the  finite  difference 
equations. 

The  detailed  development  of  the  equations  in  the  text  clarifies  the 
inherent  limitations  of  a  zero-D  model  for  an  imploding  plasma.  Besides  the 
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oversimplification  of  an  average  value  for  the  plasma  density  and  temperature, 
the  bounce  phase  of  an  imploding  annular  shell  is  poorly  handled.  Pressure 
forces  in  the  plasma  between  the  shell  and  the  axis  are  always  ignored  in  a 
zero- D  code,  but  do  become  significant  at  the  time  of  bounce.  Within  the 
framework  of  the  ZPIMP  code  this  problem  could  readily  be  corrected  by  solving 
the  dynamics  over  the  inner  zone,  as  well  as  the  plasma  shell.  Since  the  mass 
in  the  inner  zone  is  much  less  than  that  in  the  shell  zone,  one  could  still 
restrict  the  time  consuming  radiation  transport  to  the  main  plasma  shell.  For 
the  existing  ZPIMP  code,  the  results  for  general  trends  at  the  collapse  time 
are  consistent  with  several  runs  made  with  multi-zone  codes.  The  ZPIMP 
simulation  is  also  limited  to  rapid  collapse  times  since  the  gradient  scale 
length  of  the  magnetic  field  was  assumed  to  be  small  compared  to  the  plasma 
radius.  This  condition  allowed  a  simple  solution  to  the  induction  equation 
for  annular  shells.  Future  improvements  to  the  hydrodynamics  and  magnetic 
diffusion  in  ZPIMP  are  under  development. 

The  zero-D  ZPIMP  code  has  been  applied  to  the  implosions  of  an  aluminum 
vapor  in  annular  shells  and  filled  cylinders.  The  parameters  of  the  circuit 
were  chosen  to  match  the  NRL  GAMBLE  II  device.  The  effects  of  various  mass 
loadings,  initial  outer  radii,  and  shell  widths  on  the  gross  properties  of  the 
implosion  were  discussed.  In  particular,  results  for  the  implosion  time  and 
radius,  maximum  charge  state  and  temperature,  and  peak  radiated  power  have 
been  summarized  in  the  figures.  In  a  future  report  we  will  analyze  the 
radiative  yields,  spectra,  and  pulse  durations  as  a  function  of  the  mass 
loading  and  initial  configuration  of  the  aluminum  vapor.  This  large  parameter 
space  can  only  be  effectively  studied  with  a  versatile  and  rapid  zero-D  code 
such  as  ZPIMP.  From  a  comparison  of  the  results  we  intend  to  determine  the 
optimum  conditions  for  aluminum  K-  and  L-shell  radiation. 
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Uati/l«ngtli(aUro|r  aaVcta) 


1.  Minimum  outer  radius  for  the  thin  shell  A1  plasma  implosions. 
The  solid  line  represents  a  shell  thickness  AR=0.1  cm,  the  long 
dashed  line  represents  AR=0.2  cm,  and  the  short-dashed  line 
represents  AR=0.4  cm. 


2.  Time  of  maximum  implosion  for  the  thin  shell  cases.  Same 
notation  as  Figure  1. 


Maii/Leo|th(micro|f  Mn/cm) 


Fig  3. 


Maximum  electron  temperature  for  the  thin  shell  cases 
notation  as  Figure  2. 


Vat«/Len(lh(iDlcro|r  en^os) 


Fig  4.  Maximum  average  ionic  charge  for  the  thin  shell  cases 
notation  as  Figure  1. 


30  31  32  33  34 

Mass/Leng 

Fig  5.  Maximum  total  radiated  powe 
notation  as  Figure  1. 


Cotilinuum  Radiation  Below  K  Edge 


Fig  7.  Maximum  continuum  radiation  below  (a)  and  above  (b)  the  K  edge 
for  the  thin  shell  cases.  Same  notation  as  Figure  1. 


taltUI  lt>iilut(ciB) 


Fig  13. Maximum  line  radiation  below  (a)  and  above  (b)  the  K  edge  for 
the  filled  cylinder  cases.  Same  notation  as  Figure  8. 
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